On the relativistic Lattice Boltzmann method for quark-gluon plasma simulations 
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In this paper, we investigate the recently developed lattice Boltzmann model for relativistic hydro- 
dynamics. To this purpose, we perform simulations of shock waves in quark-gluon plasma in the low 
and high viscosities regime, using three different computational models, the relativistic lattice Boltz- 
mann (RLB), the Boltzmann Approach Multi-Parton Scattering (BAMPS), and the viscous sharp 
and smooth transport algorithm (vSHASTA). From the results, we conclude that the RLB model 
departs from BAMPS in the case of high speeds and high temperature(viscosities), the departure 
being due to the fact that the RLB is based on a quadratic approximation of the Maxwell- Jiittner 
distribution, which is only valid for sufficiently low temperature and velocity. Furthermore, we have 
investigated the influence of the lattice speed on the results, and shown that inclusion of quadratic 
terms in the equilibrium distribution improves the stability of the method within its domain of 
applicability. Finally, we assess the viability of the RLB model in the various parameter regimes 
relevant to ultra-relativistic fluid dynamics. 

PACS numbers: 47.11.-j, 12.38.Mh, 47.75. +f 
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I. INTRODUCTION 

Recent experiments on heavy-ion (Au-Au) collisions 
with ultra-relativistic energies at the Relativistic Heavy 
Ion Collidcr(RHIC) at Brookhaven National Laboratory, 
have revealed a new state of matter, the quark-gluon 
plasma, whereby hadronic matter undergoes a deconfin- 
ing transition which liberates quarks in the form of a gas 
of quasi-free particles [l[ ■ The quark-gluon plasma, which 
is credited for dominating the primordial state of the Uni- 
verse in its earliest 1 — 10 microseconds, shows very inter- 
esting properties, such as near-perfect fluid-like behavior, 
characterized by ultralow dynamic shear viscosity and 
associated onset of shock wave propagation 0, Q • The 
study of such shock waves plays a major role in the char- 
acterization of the quark-gluon plasma, since they carry 
information both on its equilibrium (equation of state) 
and non-equilibrium (transport coefficients) properties. 

In the last years, several numerical tools have been 
used for the computational investigation of quark-gluon 
plasmas, such as, for instance: the Boltzmann approach 
of multiparton scattering (BAMPS) [||, which solves the 
full Boltzmann equation, p tl d ll f{x,p) = C(x,p), with p^ 
the microscopic 4-momentum vector, f(x,p) the single 
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particle distribution function, and C(x,p) the binary col- 
lision term; and the viscous sharp and smooth transport 
algorithm (vSHASTA) , which is based on a hydrody- 
namic description. 

Recently, the relativistic lattice Boltzmann (RLB) 
model [f| was introduced. This new method is a rela- 
tivistic extension of the classical lattice Boltzmann (LB) 
equation^, whereby a minimal form of the Boltzmann 
equation is discretized on a lattice and the collision 
term is applied in a single relaxation time, the so-called 
Bhatnagar-Gross-Krook approximation, Ref. Q). In the 
above, "minimal form", implies that the particle veloci- 
ties are constrained to take only a very limited set of dis- 
crete values, typically of order 10 and 20 in two and three 
dimensions, respectively (sec FigJIJ, which represents an 
enormous simplification as compared to the true Boltz- 
mann equation. The key is that the proper selection of 
this handful of discrete speeds is sufficient to compute 
exactly the low-order kinetic moments which character- 
ize the hydrodynamic regime. Crucial to the success of 
this procedure is that the system be only weakly out 
of equilibrium, so that the particle distribution function 
remains close to a local Maxwellian, whose low-order ki- 
netic moments can be computed exactly with just a few 
suitably chosen discrete velocities (the nodes of Gaussian 
integration). 

The LB method shows many advantages, as for ex- 
ample, the relatively easy implementation of the simula- 
tions of fluids in complex geometries, excellent suitability 
to parallel implementation, and high flexibility towards 
the inclusion of additional physics, besides sheer hydro- 
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dynamics ^ can thus be expected that the RLB 

model would carry many of these assets over to the rela- 
tivistic context. To date, it has been shown to compute 
shock-wave formation and propagation in quark-gluon 
plasmas at a fraction of the cost of relativistic hydro- 
dynamic codes jf| [H[ . 

In order to characterize and improve the RLB model, 
further studies on its capabilities, limitations and compu- 
tational performance are in order, which is precisely the 
focus of this paper. To this purpose, it is worth reminding 
that the actual RLB model deals with weakly relativis- 
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tic fluids, characterized by 1 < 7 < 2 and £ = > 1, 

where 7 = (1 — u 2 /c 2 ) -1 / 2 is the Lorcntz factor associated 
with a fluid with speed u and T is the fluid temperature. 
However, since 7 = 2 corresponds to (3 — u/c ~ 0.85, and 
C = 1 associates with ultra-high temperatures, typically 
above 10 13 Kelvin degrees, the weakly relativistic regime 
embraces nonetheless a substantial number of interesting 
applications, including the quark-gluon plasma generated 
by recent experiments on heavy-ions and hadron jets (l5l - 
as well as astrophysical flows, such as interstellar gas 
and supernova remnants [22I [23| . 

An extension of the original RLB model, capable 
of dealing with non-Minkowskian geometries and ultra- 
relativistic fluids, has been recently developed [24j. How- 
ever, since this method is comparatively more elaborate 
than the original RLB, in the sequel we shall confine our 
attention to the latter. 

This work is organized as follows: first, in section |H] 
we give a short introduction to the RLB model, and 
show some simulations to validate our implementation. 
Subsequently, in section IIII| we extend the equilibrium 
distribution function for the particle density, to obtain 
the correct density profile in the moderately relativistic 
regime 7^2. Finally, in section IIV1 we analyze the 
effect of the lattice speed on the results, provide com- 
parisons in the moderately relativistic and high viscosity 
regimes with the previous mentioned methods, BAMPS 
and vSHASTA, and assess the viability of this scheme for 
ultra-relativistic fluid dynamics. 



II. RELATIVISTIC LATTICE BOLTZMANN 
MODEL 

In this section, we describe the relativistic lattice 
Boltzmann, which is the basis of this work, develop some 
improvements and test the new implementation against 
known results for the Riemann problem. The Riemann 
problem has piecewise constant initial conditions and a 
single discontinuity, and is commonly used to test nu- 
merical methods that model systems with conservation 
laws. 




FIG. 1. The D3Q19 (19 speeds in 3 dimensions) cell configu- 
ration. Every arrow corresponds to the scaled velocity vector 
8tBi, and the points correspond to the discretized space coor- 
dinates x. 

A. Conservation laws for relativistic fluid dynamics 

The conservation laws for relativistic fluid dynamics, 
which our RLB model is based upon, are: the conser- 
vation of particle density (in a relativistic regime mass 
is not invariant), and the energy momentum conserva- 
tion, which has to be treated separately in contrast to 
the classical LB. 

The conservation of the particle number density n is 
given by: 

d t (717) + d a (njUa) = , 

where u denotes the macroscopic fluid velocity. Here, 
Latin indices denote the three dimensional space coor- 
dinates. Note that in the following the speed of light c 
and the Boltzmann constant k are taken to be unity for 
notational simplicity. As a result, j3 = u and Q = m/T. 
Also we use the Einstein summation, i.e. repeated in- 
dices are summed upon. The hydrodynamic equations 
for the energy momentum conservation read as follows: 

d t ((e + P) 7 2 - P) + d a ((e + P) l 2 u a ) 

+ dtir 00 + d a n ab = , 
d t ((e + P) 7 2 u ab ) + d b P + d a ((e + P) i 2 u a u b ) 
+ d t TT 0b + d b ir ab = , 

where e is the energy density, P is the hydrostatic pres- 
sure, and ir ab are the components of the dissipation 
stress-energy tensor. Here the index denotes the time 
component. 



B. Implementation of the RLB scheme 

We implement the relativistic lattice Boltzmann model 
along the lines proposed in Rcfs.0, [T^ . This model is 
based on two density distribution functions, fi and gi, 
which are associated to a D3Q19 lattice Q (see Fig. [TJ, 
with 19 discrete velocities q in d — 3 dimensions. 

These velocities take just the values 0, ±q, where 
q = || is the lattice speed. Each discrete velocity c; 
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is associated with a corresponding discrete distribution 
function, fi(x,t) = f(x,t;pi) for the fluid density, and 
gi{x, t) = g(x, t; pi), with pi = m^iCi, for the fluid energy- 
momentum. 

The distributions evolve according to the following dis- 
crete BGK Boltzmann equations [8j|: 



fi (x + m,t + St) - fi (x, t) = --(fi- .O , (2a) 

T 



i 

i x i 



(4) 




4 £ 



(5) 
(6) 

(7) 



+ c^M + St) - gi(x,t) = - — (gi - g° q ) , (2b) 

T 



where r represents the local relaxation time, f^ q and g^ q 
are the equilibrium distribution functions. 

The above equations describe a two-step lattice dy- 
namics. The left-hand side encodes the free-streaming 
of the distributions along the characteristics defined by 
the discrete velocities Si. Note that since velocities are 
constant in space and time, this term is an exact lattice 
transcription of the free-streaming term in the continuum 
Boltzmann equation. This leads to major benefits for 
the computational performance of the model because, at 
variance with hydrodynamic formulations, the informa- 
tion always travels along straight lines rather than along 
space-time changing material fluid lines. 

The right-hand side describes particle collisions in the 
form of a relaxation towards local equilibria, on a time 
scale r. This is the lattice analogue of the relativistic 
Marie model p5l |. 

In order for the relativistic LB equations (Eq. (|5J)) 
to correctly reproduce relativistic hydrodynamics in the 
continuum limit, the lattice equilibria have to be designed 
in compliance with the basic number-energy-momentum 
conservation laws. As shown in Ref. [l4j], this can be ac- 
complished through an algebraic moment-matching pro- 
cedure, leading to the following expressions: 



fi= Win^i 1 + 3 



gl q = SwoP^ 2 4 



2 + cf Jfil' 
7 2 cf c 2 



,4 1 „2 

I H 



g% = 3 Wl PY { + + - 2 



(3) 

The above probability distribution functions recover 
the macroscopic values, with the ultra-relativistic state 
equation e = 3P, provided the following identifications 
are made: 



In the above, the discrete weights take the value wo = 
1/3, for | cb | = 0, Wi = 1/18 for |cj| = q, and w; = 1/36 
for |cj| = V2ci. Based on these expressions, it is readily 
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checked that J2i w i = 1 an d c l = Ej w i c ia = "3"' where 
a = x,y,z. The latter defines the sound speed, c s , and 
shows that, by stipulating q — c, i.e. the lattice speed 
equal to the light speed, the RLB supports the relativistic 
ideal equation of state c s — ^ ■ 

It is worth emphasizing that, unlike local equilib- 
ria in continuum momentum space, lattice equilibria 
arc not unconditionally positive for any value of the 
fluid velocity. This is because continuum equilib- 
ria, both non-relativistic (Maxwcll-Boltzmann) and rel- 
ativistic (Maxwell- Juettner), are irrational functions of 
both the microscopic and hydrodynamic velocity (four- 
momentum). This is no accident, but rather the result 
of the local equilibria following from an entropy mini- 
mization principle, with the entropy additivity impos- 
ing a non-rational (exponential) functional dependence 
on the collision invariants [26| . Hence, an exact tran- 
scription of such equilibria would require an infinite se- 
ries in powers of u/c s . However, since the generic n-th 
term of such an expansion involves n-th order tensors of 
the form Ei CjCj . . . c*;, it is clear that much larger sym- 
metry groups, i.e. much more discrete velocities, would 
be needed at each step of the expansion. This would 
rapidly lead to an unmanageable complexity. It is quite 
fortunate that hydrodynamics can be reproduced by en- 
suring the correct symmetry of just fourth-order tensors, 
so that, relatively simple lattices like the D3Q19, can ac- 
complish the task. Failing such fortunate circumstance, 
no LB would ever exist. 

To model shock waves in viscous quark-gluon matter, 
this scheme has to recover a special viscosity-entropy den- 
sity ratio -j , where the entropy density s is approximated 

by s = 4n~n In A, with A = ^_ and n e i = (d G = 16 

is the degeneration for gluons and T is the temperature) 
1- 

The dynamic viscosity can be computed from the re- 
laxation time according to the following expression: 
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FIG. 2. Pressure and velocity profiles for different ^ ratios, 
at t = 3.2^. The profiles show excellent agreement with the 
results in Ref. [h| . 



The term cf(r— 4£) is standard from conventional Lat- 
tice Boltzmann theory. However, in contrast to the clas- 
sical Lattice Boltzmann method, this term is pre-factored 
by jP/c 2 (c = 1 in our units) rather than by fluid density 
P- 

For the purpose of this work, r is computed with the 
initial P and 7, so as to obtain the desired - ratio. Note 
that the initial P and 7, in general for the Ricmann prob- 
lem, are functions of the coordinates, so they can lead to 
a spatially dependent relaxation time r. To avoid this, 
we have used their spatial averages. 




FIG. 3. Comparison of the particle density using both equi- 
librium functions, Eqs. Q and (©, at t = 3.2 ^i, with the 
initial conditions given in sect ion [IV A I - = 0.005, and q = 1. 
As one can appreciate, the density profiles are basically the 
same. 
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FIG. 4. Comparison of the particle density using both equilib- 
rium functions, at t = 3.2^, with the initial conditions given 
in section HVBl 2 = 0.001, and a = 1. The result obtained 
with the old (linear) equilibrium function breaks down in the 
region z £ [—1,2], whereas for the case of the new (quadratic) 
equilibrium function, the model reproduces the correct profile 
in the whole domain. 



numerical units, to 2.495 x 10~ 7 and 1.023 x 10~ 7 , re- 
spectively. The initial temperature is constant over the 
whole domain having the value 350AfeV* (in numerical 
units 0.0314), corresponding to £ ~ 3. 

The initial particle density is computed through the 
relation n = The pressure and velocity profiles at 
3.2^ for different 2 ra tios between 0.001 and 0.1, are 

c s ' 

shown in Fig. PZl from which excellent agreement with the 
results in Ref. [141] is readily appreciated. Each simulation 
took around one second on a Intel core i5 of 2.3GHz. 



C. Numerical validation 

In order to test our numerical scheme, we carry out 
some simulations of shock waves in quark-gluon plasma [|J 
in one dimension, and compare the results with the ex- 
isting literature, Ref. 1J]. For this simulation, a lat- 
tice with 1 x 1 x 800 cells, open boundaries along the 
mainstream z-direction and periodic boundaries along 
the cross-flow a;-and y-dircctions are used. As a result, 
we set Sx = 0.008/m and St = 0.008^;. The initial con- 
ditions for the pressure arc P(z < 0) = Po = 5.43y^£ 

and P(z > 0) = Pi = 2.22^£. This corresponds, in this becomes negative for counter-streaming populations 



III. 



IMPROVING THE EQUILIBRIUM 
DISTRIBUTION 



Inspection of the particle density equilibrium dis- 
tribution function, f° q =Win"/ ^1 + 3^^-^ , shows that 
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when (ci ■ u) < —-3-. This is non-physical and leads to 
cither unreliable results (see Fig. |4|) or numerical insta- 
bility. 

In previous works [(| [T3 [. this problem was circum- 
vented by increasing the lattice speed q and reducing 
the time step accordingly, so that the light-cone condi- 
tion ciSt = 5x, remained fulfilled. In the low- viscosity 
regime this was found to stabilize the numerics without 
affecting the physics to any appreciable extent. However, 
as shown in the sequel, in at higher viscosities, this is no 
longer the case. 

A detailed analysis of the effect of the lattice speed on 
the physical results is presented in section IWl However, 
prior to this, we first introduce a new particle density 
equilibrium function, which allows larger fluid velocities 
at a given value c; = 1, without violating the positivity 
condition, /? q > 0. 

This new equilibrium distribution reads as follows: 



9(c; 
2 



(9) 

This expression includes a new quadratic term in the 
macroscopic fluid velocity. This new term is standard 
in classical lattice Boltzmann theory, and results form 
second order expansion in the fluid velocity of a locally 

TO(«-lO a 

shifted Maxwellian distribution ~ e . 

To check that the new equilibrium function is valid 
and reproduces the same results as the old one, we carry 
out several simulations in the weakly relativistic regime, 
which each of them required around one second on an In- 
tel core i5 of 2.3GHz. For illustration purposes, we show 
just one example, where we only consider the particle 
distribution function, because the improved equilibrium 
function has no direct influence neither on the pressure 
nor on the velocities, see Eq. (|4|. 

From Fig. [3l which shows that density profile at rj/s = 
5 10~ 3 , we observe that both equilibrium distribution 
functions give basically the same result. However, from 
Fig. |4l which refers to rj/s = 10 -3 , we can see that the 
computation of the particle density using the old equilib- 
rium distribution function, in the moderately relativistic 
regime, is breaking down in the region z € [— 1,2] /to, 
and leads to very large fluctuations, with both positive 
and negative values. For visualization purposes, just the 
physically meaningful region is plotted. Note that out- 
side of the region z G [—1,2] /to both equilibrium func- 
tions reproduce the same result. However, in the course 
of the evolution, those fluctuations are found to propa- 
gate over the entire simulation domain. 



IV. INFLUENCE OF THE LATTICE SPEED 

Having illustrated the stabilization effects of the new 
equilibrium distribution, all following numerical exper- 
iments are performed with this distribution. Next, we 
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FIG. 5. Weakly relativistic case, 2 — 0.001. Pressure and 
velocity profiles for different lattice speeds at t = 3.2^. In 
this case, the choice of the lattice speed does not affect the 
result of the simulation. 



study the influence of the lattice speed, cj. For this pur- 
pose, we simulate various shock waves in quark-gluon 
plasma with different viscosity-entropy density ratios and 
different lattice speeds in the weak (1 < 7 < 2) and 
mildly (7 ~ 2) relativistic regimes. As noted before, 
by increasing the lattice speed, we need to decrease the 
time step, in such a way that x + Cidt corresponds to a 
grid point in the lattice. Basically, in numerical units, 
q = 10 implies dt = 0.1, so that the product is always 
1 or depending of the velocity vector Cj. Clearly, 
smaller time-steps imply a correspondingly larger num- 
ber of time steps, hence more simulation time. 



A. Weakly relativistic regime, 7< 2 

For the weakly relativistic simulations, we use a lattice 
with 1 x 1 x 800 cells, open boundaries in z-direction 
and periodic boundaries in ;c-and y-direction. We set 
Sx = 0.008/to and St = The initial conditions 

for the pressure arc P(z < 0) = Po = 5-43y^ and P(z > 
0) = P\ = 2.22y^. This corresponds in numerical units 
to 2.495 x 10~ 7 and 1.023 x 10~ 7 , respectively. The initial 
temperature is constant over the whole domain, namely 
400Afey(in numerical units 0.036), corresponding to £ ~ 
2.5. The initial particle density is calculated with the 
equation of state, n = ^. 

A snapshot of the pressure and velocity profiles at 
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FIG. 6. Weakly relativistic case, 2 — 0.005. Pressure and 
velocity profiles for different lattice speeds and , at t = 3.2^. 
Here, some differences are visible. 
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FIG. 7. Weakly relativistic , 2 — 0.01. Pressure and velocity 
profiles for different lattice speeds at t = 3.2^p. Here we 
observe a small effect of the lattice speed. 



t = 3.2^-, for different lattice speeds(c/ = 1, cj = 10, 
and ci = 100) and different viscosity-entropy ratios (be- 
tween 0.001 and 0.01), is shown in Figs. [5j[6]and[7l From 
these figures, we can observe that , as expected based on 
previous experience, there is no significant effect of the 
lattice speed on the results. However, there is a signif- 
icant difference in the computational performance, the 
simulations took around 1,9, and 90 seconds for q = 1.0, 
q = 10.0 and q = 100.0, respectively. 



B. Moderately relativistic regime, 7 ~ 2 

For the moderately relativistic regime, we perform the 
simulations on a lattice with 1 x 1 x 1600 cells, in or- 
der to cover a larger computational domain. The same 
boundaries are used as before. We set Sx = 0.008/m and 



St 



_ 0.008 fm 



The initial conditions for the pressure are 



P Q = 5.43f^ anc i p = 0.339^. In nume rical units, 

u fm 6 1 fm 6 ' 

they correspond to 2.495 x 10~ 7 and 1.557 x 10~ 8 , re- 
spectively. The initial temperature is now different for 
z > namely T x = 200MeV (in numerical units 0.018), 
corresponding to £ ~ 5. 

For z < the initial temperature is the same as in 
the previous simulation, To = 400AfeV^ (in numerical 
units 0.036), and the initial particle density is computed 
in the same way as before. The velocity and pressure 
profiles at t = 3.2^, for different lattice speeds(q = 1, 
q = 10, and c/ = 100) and different viscosity-entropy 
density ratios (0.001 and 0.05), are shown in Figs. [8] 
and[5] Again, no significant influence of the lattice speed 
is observed. 

For 2 = 0.01, we perform the simulation using a lat- 
tice with 1 x 1 x 800 cells, so that Sx = 0.016/to and 
5t = °'"; 16 ^7? ■ In this case, the numerical units of the 
pressure are P = 1-996 x 10~ 6 and P 1 = 1.2456 x 10~ 7 , 
respectively. The numerical values of the temperatures 
and particle density are set up in the same way as be- 
fore. The results of this simulation are shown in Fig. 
[T7J1 Here, we see that the lattice speed drastically af- 
fects the results. The implemented simulation in this 
section spanned around 1, 4, and 44 seconds for the cases 
C[ = 1.0, Cj = 10.0, and c; = 100.0, respectively. 

Summarizing, we conclude that the physical results in 
the moderately relativistic regime are affected by the in- 
creasing value of the lattice speed beyond c; = 1.0. Also, 
we see that this effect is more pronounced with increas- 
ing 2 . Furthermore, we observe that the results of the 
simulation with high lattice speeds, for a given ratio 
are similar to the ones obtained with c/ = 1 for a higher 
viscosity-entropy density ratio. To understand this ef- 
fect, a deeper theoretical investigation of the basic RLB 
scheme is required. According to the Chapman-Enskog 
expansion, the relation for 77 results from a Taylor expan- 
sion of the discrete lattice Boltzmann equation to second 
order in space and time, combined with a first order per- 
turbation in the Knudsen number K n ~ c s t/8x = ^r/St, 



7 





FIG. 8. Moderately relativistic case, 2 — 0.001. Pressure 
and velocity profiles for different lattice speeds at t = 3.2 — . 
Very small differences appear only at the highest /3. 



FIG. 9. Moderately relativistic case, 2 = 0.005. Pressure 
and velocity profiles for different lattice speeds at t = 3.2 — . 
Here we see some differences for different lattice speeds. 



and quadratic truncation of the local equilibria in the 
Mach number Ma = u/c s = -s/3/3- Since, as we have 
observed before, raising c\ implies lowering St, it is clear 
that simulations with q > 1 at a given r, imply larger 
values of the Knudsen number, thus pushing RLB poten- 
tially outside of the domain of validity of the Chapman- 
Enskog asymptotics. Similar problems are well known to 
occur in the simulation of complex non-relativistic fluids 
with sharp interfaces (27| . 



V. COMPARISON WITH BAMPS AND 
VSHASTA FOR HIGH VALUES OF 2 

3 

In this section, we perform simulations to compare the 
RLB model with BAMPS and vSHASTA gj for the case 
of high viscosities and high speeds. The initial conditions 
are the same as in section HVBI We use q = 1.0 and a 
lattice with 1 x 1 x 1600 cells, so that the numerical units 
stay the same. Each simulation took around one second 
on an Intel core i5 of 2. 3 GHz. The results of the simula- 
tion at t — 3.2^ for two different viscosities are shown 
in Figs. HI] and [fe] Due to the fact that BAMPS is the 
solution of the Boltzmann equation with the complete 
collision term, while the RLB model is a near-equilibrium 
approximation, and vSHASTA is a second order scheme 
for modeling relativistic fluid dynamics, it is very inter- 
esting to compare the three techniques with each other. 

For 2 = 0.1 the results from BAMPS and vSHASTA 
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FIG. 10. Moderately relativistic case, 2 = 0.01. Pressure 
and velocity profiles for different lattice speeds at t = 3.2 — . 
Here we observe, that especially for the highest /3, the results 
differ substantially. 
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FIG. 11. Pressure and velocity profiles for the different mod- 
els with 2 = o.l. BAMPS and vSHASTA show only tiny 
departures from each other, but differ significantly from RLB. 



FIG. 12. Pressure and velocity profiles for the different models 
with 2 = 0.5. vSHASTA has stability problems and discon- 
tinuities, while the RLB diverges from the BAMPS solution 



even more. 



are very close, and the RLB presents some small devia- 
tions. On the other hand, for 2 = 0.5, BAMPS still re- 
produces stable results, whereas vSHASTA suffers some 
stability problems and discontinuities, and the deviations 
of the RLB become more pronounced. To understand 
these deviations in high viscosity regime, we start first 
by revisiting the meaning of the viscosity-entropy den- 
sity ratio within the RLB model. The viscosity-entropy 
density ratio can be written as 



(10) 



and using the relation T = ^, 



V 
s 




(11) 



3 + | In 



This expression shows that, for high temperatures, there 
is a nearly linear dependence of 2 on the temperature, 
which means that when we simulate a high 2 ratio, this 
corresponds, in terms of the Boltzmann equation, to a 
high equilibrium temperature. 

Note that the results of the RLB simulations present 
a discontinuity at z — fm in Figs. [TT] and Q21 This 



discontinuity appears as a consequence of the initial con- 
dition for the pressure (solid line in Fig. [T0|) . which is 
set according to the Riemann problem. Therefore, we 
can conclude that the RLB, in the case of high viscosity- 
entropy density ratios, is not able to solve correctly the 
Riemann problem and maintains the initial discontinuity 
during the whole simulation. 
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FIG. 13. Maxwell- Jiittner distribution for different velocities 
at very low temperatures, mc 2 /kT = 200, the dashed lines 
denote the parabolic approximation to the respective MJ dis- 
tribution, denoted by solid lines. Only the positive velocities 
are shown since the distribution is symmetric. 
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The RLB model approximates the probability equilib- 
rium distribution functions by a quadratic expression, i..e 
a parabolic approximation. At sufficiently low tempera- 
ture, this approximation is accurate, but as temperature 
is increased, hence higher values accuracy is rapidly 
lost. To appreciate this effect, it is instructive to graph- 
ically inspect the Maxwell- Juttner distribution (MJ)[28j| 



fMj(x,p,t) 



(-*?*). 



which depends on the 



z CX P ^ T 

macroscopic four- velocity U a = j(u) (1, u), and the mi- 
croscopic four-momentum p a = ^{v)m (1, —v), where m 
is the rest mass and I; the microscopic velocity and Z 
a normalization factor that depends on the temperature 
and macroscopic velocity. 

From Figs. Q21 HH and [THl we observe that the MJ 
distribution becomes broad and very asymmetric for 
high velocities and high temperatures, and therefore a 
parabolic approximation is no longer valid. In fact, at 
C < 1, the parabolic approximation should be replaced by 
a bimodal parabolic one, which is certainly feasible, but 
beyond the current RLB formulation. This explains why 
RLB has problems to reproduce results in this regime. 
However, for the case of low temperature and relatively 
high velocities, the RLB model can still model the proper 
fluid dynamics since its equilibrium distribution func- 
tions (parabolic approximation) remain very close to the 
MJ distribution. For high temperatures, it is only possi- 
ble to study weakly relativistic systems (see Fig. [T5)) . 
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FIG. 14. Maxwell- Juttner distribution for different velocities 
at low temperatures, mc 2 /kT = 20, the dashed lines repre- 
sent the parabolic approximation to the respective MJ distri- 
bution, denoted by solid lines. 

As a general conclusion, we expect that the RLB ap- 
proach can model moderately-relativistic fluid dynamics 
(\u\ « 0.9) only in the low temperature regime, £ > 1. 
For illustration purposes, we show a qualitative sketch of 
the regimes characterized by the parameters £ and 7 in 
Fig. [TO] 
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VI. CONCLUSION 

Summarizing, in this paper we have determined the 
capabilities and limitations of the RLB model for the 
simulation of relativistic flows in the various regimes as- 
sociated to the flow speed and temperature. For this pur- 
pose, we have performed extensive simulations of shock 
waves in quark-gluon plasma at weakly and moderately 
relativistic regimes and low and high viscosities. We 
have introduced a higher order equilibrium distribution 
function, which is shown to enhance the stability of the 
original RLB model, and permits to compute correctly 
the particle density profile in the moderately relativis- 
tic regime, with no need of enhancing the lattice speed, 
hence no need of reducing the time-step, as in the orig- 
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inal RLB model. The result is that RLB can compute 
well-resolved weakly and moderately one-dimensional rel- 
ativistic shock wave propagation in less than a minute, 
CPU time on an Intel core i5 of 2.3Ghz. To the best of 
our knowledge, this is significantly faster than any hy- 
drodynamic code, let alone the full BAMPS solution. 

Furthermore, we have performed a numerical investi- 
gation of the effect of the lattice speed on the physical 
results, and concluded that the results can differ, due to 
higher order terms in the dissipation tensor which are 
not included in the Chapman-Enskog analysis. As a con- 
sequence, we observed that the results of the simulation 
with high lattice speeds, for a given ratio -, are similar 
to the ones obtained with c; = 1 for a higher viscosity- 
entropy density ratio. 

In this study, we also showed that the choice of high 
viscosity-entropy density ratios at high speeds, affects the 
accuracy of the results. In addition, by direct comparison 



of the current RLB model with BAMPS and vSHASTA 
in a moderately rclativistic and highly viscous regime, 
we have shown that the parabolic approximation of the 
Maxwell- Juettner distribution poses significant restric- 
tions to the viability of RLB for strongly relativistic, high 
temperature, fluids. 

Finally, based on the study of the approximation of the 
MJ distribution with parabolic equilibria, we are led to 
predict, that the current RLB model would properly work 
at relatively low temperatures, £ > 1, and high velocities 
(f3 ~ 0.9). Fortunately, such regimes are by no means 
devoid of interesting physical applications. Extension to 
more general relativistic flows requires further develop- 
ments, such as the introduction of higher order lattices, 
with higher order equilibria, possibly equipped with rel- 
ativistic H-theorems (entropic LB methods). This offers 
a very interesting object of future research in the field. 
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